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ABSTRACT 

We perform a dynamical analysis of the recently updated set of the radial velocity (RV) measurements of 
the HD 160691 (/j Arae). The purely kinematic, 2-Keplerian model of the measurements leads to the best-fit 
solution in which the eccentricity of the outer planet is about 0.7 and its semi-major axis is about 4 AU. The 
parameters of the inner planet are well determined. The eccentricity is about 0.3 and the semi-major axis is 
about 1.65 AU. In such circumstances, the best 2-Keplerian model leads to a catastrophically unstable config- 
uration, which is self disrupting in less than 20,000 yr. In order to derive dynamically stable configurations 
which are simultaneously consistent with the RV data, we use the so called GAMP (Genetic Algorithm with 
MEGNO Penalty) which incorporates the genetic fitting algorithms and a verification of the stability by the 
fast indicator of the dynamics (MEGNO). Using this method, we derive meaningful limits on the parameters 
of the outer planet which also provide a stable behavior of the system. It appears, that the best-fit solutions 
are located in a shallow valley of (xl)^^^, in the (flc,ec) -plane, extending over 2 AU (for the formal lo confi- 
dence interval of the best fit). We find two equally good best-fit solutions leading to the qualitatively different 
orbital configurations. One of them corresponds to the center of the 5:1 mean motion resonance (MMR) and 
the second one describes a configuration between the 6:1 and 17:2 MMRs. Generally, the HD 160691 system 
can be found in the zone confined to other low-order MMRs of the type p: 1 with p > 5. Our results support the 
classification of the /j Ara as a hierarchical planetary system, dynamically similar to other known multi-planet 
systems around HD 12661 and HD 169830. The results of the GAMP analysis on the extended data set are 
fully consistent with our previous conclusions concerning a much shorter observational window. It constitutes 
a valuable confirmation that the applied method is well suited for the analysis of the RV data series which only 
partially cover the longest orbital period. 

Subject headings: celestial mechanics, stellar dynamics — methods: numerical, N-body simulations — ^planetary 
systems — stars: individual (HD 160691) 



1. INTRODUCTION 

Precise measurements of the radial velocity (RV) of two 
stars, HD 160691 (jj Ara ) and HD 154857, from the Anglo- 
Australian survey for extrasola r planetary sy stems have been 
recently published by McCarthv et al. ( 2004). O ur particular 
attention is devoted to the HD 160691 system. JJutler et al.l 
(^001) reported a planet with a 2 yr orb ital period about the 
jU Ara . This discovery was confirmed bv lJones et aL (i2002b) 
and the observational window extended to 4 yr revealed an 
additional linear trend in the RV signal, interpreted as the sig- 
nature of a more distant planet. The data presented in the 
new paper by the Anglo- Australian Planet Search team cover 
about 6 yr This extended observational window allows to 
confirm the existence of the seco n d Jo vian planet in the sys- 
tem. Very recently, Santos et al. (2004) discovered a small, 
14 -Earth masses object in a circular orbit, at the distance of 
onl y 0.09 AU from the p Ara star and confirmed the findings 
of dMcCarthv et alJ2004l) . The Anglo- Australian team deter- 
mined the Keplerian elements of the Jovian companions: the 
orbital periods of 645 d, 3000 d and the eccentricities of 0.2, 
0.57, respectively ( McCarthv et al, 200 4). According to these 
authors, the current data set covers about 70% of the longer 
period. 

The detection of multiple, Jovian-like planets in stages from 
shorter to longer orbital periods are obviously expected from 
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the Doppler surveys and dMcCarthv et alJ2004^ announce an- 
other possible multiple system around HD 154857. It is in- 
teresting to follow the subsequent estimations of the orbital 
parameters of such system, derived on the basis o f an increas- 
ing number of observations. In the preprint by 'Jones et alj 
(|2002a), the preliminary orbital solution for /j Ara suggested 
a proximity of the two putative planets to the 2:1 mean mo- 
tion resonance (MMR). Such configuration would be the third 
occurrence of this particular low-order resonance among only 
about 15 multi-planet systems known at that time. However, 
the formal 2-Kepler solution found in the preprint appeared to 
be very unstable. In the published paper, Jones et al. ( 2002hJ) 
excluded the possibility of the 2:1 MMR due to an increased 
orbital period of the outer planet (by abou t 300 d). Our inde- 
pendent analysis of the same RV data set dGozdziewski et all 
2003) has revealed that in the specific case of /j Ara , when one 
has only a part of the orbital RV signal of the additional plane- 
tary companion, the application of 2-Kepler model of RV sig- 
nal certainly leads to artifacts. The formally best-fit configura- 
tions are self-disrupting on a short time-scale! The newest set 
of measurements leads to the 2-Keple r solution with the oute r 
planets' period increased to ~ 3000 d ( McCarthy et al. 2004^). 
This 2-Keplerian fit still describes a configuration which dis- 
integrates after only 20,000 yr. It can be easily checked by 
direct integration of the 3-body problem. In fact even a sim- 
ple analysis shows that this best fit is located in the proximity 
of the planetary collision zone. The exact collision curve can 
be determined through a\,{l +e\,) ~ 0^(1 — Cc)- This means 
that the large eccentricity of the companion c would be pos- 
sible only if an orbital dynamical mechanism protecting the 
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planets from an encounter existed. Another possibility is that 
the real o rbital elements are in fa ct different from those deter- 
mined bv' McCarthv"etal] ( l200l . 

Clearly, one needs a method of dynamical analysis of the 
orbital fits which will reproduce an observed RV signal and 
simultaneously fulfill stability criteria. In our method, we not 
onl y account for the mutual inte ractions between the plan- 
ets jLaughlin & Chambers! I200I but also eliminate the or- 
bital fits corresponding to strongly chaotic orbits, regardless 
of their goodness of fit expressed by (Xv)'''^- We treat the 
dynamical behavior in terms of the chaotic and regular (or 
weakly-chaotic) states as an observable, at the same level 
of importance as the RV measurements are. Details of the 
algorithm which will be called from hereafter the Genetic 
Algorithm with MEGNO Penalty (GAMP) are presented in 
JGozdziewski et al. 2003 ). This method seems to be particu- 
larly useful for systems undergoing strong mutual interactions 
such as the /j Ara system. According to the quoted prelimi- 
nary estimates of the semi-major axes of the companions and 
their eccentricities, the likely configurations are related to the 
zone of strong low-order mean motion resonances like 5:2, 
4:1, 5:1, 6:1, 7:1. The widths of these resonances rapidly in- 
crease with the growing eccentricity of the outer companion, 
finally leading to their overlapping and a creation of a zone 
of global instability (see section 4 in this work). In this in- 
stance, the MMRs-originating chaotic diffusion leads to vio- 
lent, macroscopic change s in the physical state of the system 
iFroeschle & Legalll999ft . In a planetary system, it means 
close encounters between the planets and/or the parent star, 
collisions or ejection of a body from the system. Following 
this reasoning, to find configurations which are stable over 
long periods of time, it is essential to eliminate the strongly 
chaotic solutions. These solutions originate in the unstable 
MMRs and lead to a discontinuity of the searched parameter 
space, in the sense of the dynamical stability. The sophisti- 
cated structure of the phase space can be understood in terms 
of the KAM theorem. Obviously, both the kinematic models 
of the Doppler signal as well as the gradient-like methods of 
the RV data fitting which require the (Xv)^''^ to be a continu- 
ous function of its variables, are not well suited to this case. 
On the contrary, the discontinuous structure of the orbital pa- 
rameter space justifies the application of non-gradient fitting 
and Monte-Carlo based techniques. 

Taking into account these very ba sic considerations, we an- 
alyzed the RV data published by Jones et al. (2002b). We 
found that the quasi-linear trend present in the reflex motion 
of HD 160691 permitted a continuum of equally good solu- 
tions in the (flc,ec) -plane. Moreover, using the GAMP algo- 
rithm and assuming that <5 AU, it was possible to reduce 
the space of stable fits, permitted by the RV data, to a num- 
ber of solutions confined to the zone of low-order mean mo- 
tion resonances (MMRs) like 3:1, 7:2, 4:1, 9:2, 5:1, 11:2, 6:1. 
Moreover, the short observational window did not permit a 
distinction between them. At present, having the access to 
the updated observations, we have an occasion to verify these 
early predictions and to test the reliability of the fits obtained 
with GAMP. Many stars accompanied by Jovian planets re- 
veal linear trends in the Doppler signal indicating the presence 
of additional companions (e.g., Butler et al. 2003). Hence, we 
have a chance to test if our method is helpful in obtaining pre- 
hminary and realistic estimates of their orbital parameters. 

2. THE NUMERICAL SETUP 



As the basic tool for searching for the best fit solutions 
to the RV data, we use the genetic a lgorithm scheme (GA) 
implemented by ICharbonneaul Jl995h in his publicly avail- 
able code named PIKAIA^. The synthetic reflex motion of 
the star is described by the 2-Keplerian RV model (as in 
iGozdziewski et al ."2003l) and the self-consistent Newtonian, 
A^-body model (Laughlin & Chambers 200 11) incorpor ating 
the MEGNO indicator penalty (.Gozdziewski et alJl20'03i) . It 
is known that the genetic scheme is not efficient for finding 
very accurate best-fit solutions but provides good starting fits 
for much faster and more precise gradient methods, like the 
Levenberg-Marquardt (LM) algorithm (Stepinski et al. 2000). 
Nevertheless, the GA makes it possible, in principle, to find 
the global minimum of (xl ) ' while the gradient methods are 
basically local. In our experiments, the GA fits were finally 
refined by a yet another very accurate non-grad ient m inimiza- 
tion scheme by Melder and Mead ( Press et al1 ll992l) . widely 
known as the simplex method. The fractional convergence 
tolerance to be achieved in the simplex code has been set in 
the range 10^^-10^^. Typically, the lower accuracy has been 
forced by time consuming GAMP tests. The application of 
the simplex refinement in the CPU expensive GAMP code is 
essential because it reduces the CPU usage dramatically (by a 
factor of ten, at least). 

For the stability analysis, we use two numerical tools which 
are helpful in resolving the regular (quasi-periodic) and ir- 
regular (chaotic) character of a teste d configuration. The al- 
ready mentioned MEGNO algorithm (Cincotta & Sim6"2000'; 
Cincotta et al. 2003; Giordano & Cincotta 2004) is already 
used in a few our papers. In this work (also to avoid a rou- 
tine approach in our stability studies), we apply an alterna- 
tive method derived on the basis of the Fast Fourier Trans- 
form (EFT). This spectral method (SM) was developed by 
Michtchenko & Ferraz-Mello (2001). Its idea is very sim- 
ple: in the quasi-periodic system, the spectral signal, (here 
/(?) = fl(f)exp(iA,(f) where a(t) and X(t) are the heliocen- 
tric canonical elements ( Laskar & Robutel 1995), the semi- 
major axis and mean longitude, respectively), has only a small 
number of leading frequencies while in the chaotic system 
the spectrum of these frequencies becomes very complicated. 
By counting the number of the frequencies (the spectral num- 
ber, SN) in the EFT spectrum of the signal that have ampli- 
tudes over a given noise level (typically, we use a threshold of 
about a few percent of the largest amplitude), one can distin- 
guish between a regular and chaotic motion. It appears that 
these two fast indicators are similarly well sensitive to the 
chaotic diffusion generated by the MMRs in the systems with 
Jupiter-like planets and the required integration time is rela- 
tively very short, typically, about lO'* periods of the outermost 
planet. The spectral method seems to have many advantages. 
Its simplicity is appealing. Under certain conditions, the SM 
is even more efficient than MEGNO because one avoids inte- 
grating complex variational equations. It provides a straight- 
forward identification of the MMRs. To summarize the dif- 
ferences between the two methods, one could say that the SM 
detects the chaotic diffusion as an "analogue" device, while 
the MEGNO (Lyapunov exponent type indicator) is more of 
a "digital" tool. Details of the comparison of these two com- 
plementary methods will be described elsewhere. 

3. SEARCHING FOR THE NEWTONIAN BEST FITS 
'' ihttp ; / /www . hao . ucar . edu/public/ research/ si /pikaia/pikaia . html| 
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First, we performed the analysis of the 2-Keplerian fits 
to the RV data of the HD 160691 system similarly to 
(|Gozdziew ski et alJ200'3h . The best fit solution has (Xv)'''^ = 
1 .644 and the rms of 4. 147 m/s (see Table 1). This fit is catas- 
trophically unstable. It is not surprising because ~ 0.72 by 
far exceeds eccentricities corresponding to the collision curve 
for the two planets. Our 2-Keplerian fit has significantly larger 
gc and the orbital period of the outer planet than found by the 
discovery team. At present, Pc — 3360 d is also longer than 
from the previous 2-Keplerian estimate (tTones et al.ll2002bt 
iGozdziewski et al]l2003i) . This conclusion and our experi- 
ments (not only from this work) seem to follow an empirical 
rule: if the measurements cover only a fragment of the longest 
orbital period, widening the observational window leads to an 
even longer period. Obviously, the current set of RV data still 
does not permit to identify a stable solution, when assum- 
ing the purely kinematic 2-Keplerian model of the RV data. 
The synthetic RV curve of the best 2-Keplerian fit in Fig.^is 
barely different from the curve corresponding to the synthetic 
curve derived from the A^-body model (more details will be 
given further). Yet, they describe completely different config- 
urations. Obviously, their (Xv)^^^ rms must be almost the 
same. This example clearly shows that the smallest (Xv)'''^ or 
rms cannot be the decisive factor for the choice of the proper 
description of the observed planetary configuration. 

In the next experiments, we used only the full A^-body 
model of the reflex motion of the parent star The first ques- 
tion is how (if at all) the self-consistent fits (for the moment — 
without the instability penalty) change the character of the 
best fit solutions. Such fits have been initialized by the GA 
and refined by the simplex. The results are illustrated in 
Fig. 12] It shows the projections of the best fits to which the 
simplex converged in every individual run (expressed by os- 
culating heliocentric elements) on the planes of the relevant 
orbital parameters. Hundreds of independent runs of the fit- 
ting code have been preformed. The osculating elements are 
given at the moment of the first observation. The best fit so- 
lution is given in Tabled as the NLl fit. Its rms is about 
4.15 m/s and {xiy^^ — 1.643, thus almost the same as for 
the 2-Keplerian model. Yet, the derived minimal mass of the 
outer planet is much smaller than in the 2-Keplerian solution. 

The collected fits make it possible to obtain a general view 
of the orbital configurations of the system permissible by the 
A^-body model. It appears that the elements of the inner planet 
are already well fixed. This is not the case for the outer 
companion. Its eccentricity, in the range of lo confidence 
interval of the best fit solution NLl, can vary over 0.1-0.8 
and the semi-major axis spreads over 2-3 AU. Also the mass 
of the outer planet cannot be precisely estimated. Surpris- 
ingly, the angular elements of both companions are relatively 
well bounded. Assuming a coplanar configuration, this re- 
sult makes it possible to reduce the space of barely deter- 
mined parameters to three dimensions: {Mc,ac,ec). This is 
quite comfortable for the studies of the system dynamics be- 
cause it can be limited to a representative {ac,ec)-p\ane pa- 
rameterized by the mass of the outer companion. The as- 
sumpti on of coplanarity may be justified by the recent re- 
sults (Santos etal. 2004). As we already mentioned, these 
authors discovered a small, Neptune-like planet around the 
HD 160691 . They also measured the projection of the ro- 
tational velocity, vsin/, of HD 160691 and found it close to 
the theoretical, maximal value. Let us suppose that the orbital 
plane of the smallest planet is almost perpendicular to the ro- 



tational axis of HD 160691 . This means that the orbit is al- 
most "edge-on". In the presence of a possibly highly inclined 
Jovian planet, the orbit of the smaller companion would be 
influenced by the Kozai resonance. For a large relative incli- 
nation, this mechanism forces a large eccentricity of the inner, 
smaller body and finally can destabilize its motion. Because 
the ec centricity of the inner planet is very small (pantos et aj 
12004ft . the two most inner orbits should not be significantly 
incUned. Assuming that the planetary system emerged from 
a coplanar proto-planetary disk, the whole system is likely 
close to the coplanar configuration. Still, the recent results by 
Thommes & Lissauer (20QJ and^d ams & Lauahlin (20031) 
indicate that a large fraction of planetary systems which inter- 
act gravitationally can be in fact significantly non-coplanar 
Unfortunately, the currently available data do not permit to 
investigate such a possibility. 

As we have akeady noticed, the NLl solution is formally 
as good as the fit obtained with the 2-Kepler model. Un- 
fortunately, there is also a similar problem with the stabil- 
ity of the relevant configuration. The evolution of the orbital 
elements (Fig. |3} clearly shows that the motion is strongly 
chaotic and, in fact, corresponds to collisional orbits. In the 
specific integration, the Lyapunov time is only 233 yr. It 
will be shown later that this fit corresponds to the border of 
5:1 MMR (Fig.[lOj. 

In order to estimate t he formal err ors of the best fit NLl, 
we followed Gozdziewski & Macieiewski (2003). The actual 
error of the determination of the RV consists of two parts: the 
measurement error and the contribution from the variability 
of the star itself. For chromospherically quiet G and K dwarfs 
the stallar jitter is about the level of the measurements eiTors. 
In this paper we adopted Ojittei — 3 m/s, fol lowing the esti- 
mates of the stellar variability of G dwarfs bv lMcCarthv et all 
(2004) (see their Fig. 1). These estimates are based on 6 yr 
monitoring of a few representative G-type stars similar to 
/J Ara . To estimate the error of the best fit NLl, we added 
to every original measurement a random Gaussian noise with 
the zero mean value and the mean dispersion o ~ 4 m/s and 
the Gaussian noise of the stellar jitter with Ojittei — 3 m/s. 
Next, the best fit was recalculated with the simplex algorithm 
using the synthetic data set. This procedure was repeated 
500() times. After finding the mean values of the best fit pa- 
rameters, we calculated their dispersions. These dispersions 
are used as the estimates of the mean errors of the osculating 
elements and shown in parentheses in Table 1. These esti- 
mates confirm our earlier findings indicating that the plane- 
tary phases can be relatively well fixed while the mostly un- 
constrained elements are the mass, the semi-major axis and 
the eccentricity of the outer planet. 

Further, to see whether the formal RV eiTors increased by 
a contribution from the stellar jitter influence the overall dis- 
tribution of the best fit elements, we repeated the GA search 
using the original data points and their mean errors changed to 

^ ^ y'c?^ + C^j^itter- parameters of the best fit found in this 
search are given in Table 1 as the NL2 solution and the distri- 
bution of the best fit elements is shown in Fig. |3 and Fig.ls] 
Obviously, the larger mean errors do not lead to a significant 
modification of the overall distribution of th e best fits. Let us 
note that the projections shown in Fig. I2I4I and Isl serve as an 
alternative, realistic estimation of the fit errors. The eiTors ob- 
tained in this way are even larger than the previously derived. 

The experiments utilizing the A^-body model of the RV mea- 
surements confirm that the existence of the second planetary 
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companion. Nevertheless, all the best-fit initial conditions still 
lead to self-disrupting configurations and, additionally, are not 
well constrained. Recalling the discussion in the introduction 
and th e results of dynamical analysis in Gozdziewski et al. 
(|2003), we already have a rough picture of the phase space in 
the relevant zone of (0^,6^). This zone is confined to strong, 
unstable low-order MMRs. Due to the preferably large initial 
Cc, the chances of finding stable configurations by a straight- 
forward Keplerian or Newtonian fit to the RV data are poor. 

Thus, to find the best-fit orbital solution in a zone of sta- 
bility, we have made extensive calculations by applying the 
GAMR The code has been restarted hundreds of times. The 
control parameters of the code have been changed and ad- 
justed to give an additional degree of randomness in the GA 
search. The results are illustrated in Fig.|6l As in the previ- 
ous cases, this figure shows a projection of the found best-fit 
solutions onto different planes of osculating elements. Only 
those fits are shown which have (Xv)'''^ < 1.668. This limit 
corresponds to the formal la confidence interval of the best 
fit solution (GMl, see Table 1) having {xlY^^ ^ 1.64701 
and rms 1.147 m/s. As we can expect, the stable solutions 
have er < 0.5. This stro ngly confirms our prediction from 
iGozdziews ki et al.l ( l2003h . 

The reader may notice that the (Xv)'''^ for GMl (Table 1) 
is given with an unusually high precision of 5 decimal places. 
This is due to another solution (GM2, see Table 1) having 
almost the same {xiy^^ while its Uc is shifted by 0.4 AU. 
The orbital behavior of the corresponding configurations (see 
Fig. and Fig.|8j clearly shows that they belong to qualita- 
tively different dynamical regimes. The GMl fit corresponds 
to a very stable, regular behavior which is confirmed by the 
MEGNO rapidly converging over 1 Myr Surprisingly, the 
GM2 fit describes a system whose stability is marginal. This 
system is trapped into the 5:1 MMR. This is illustrated in 
panel Fig. 0; showing one of the critical arguments of the 
resonance, Osa = SX^ — — (S^ — 303c. 

4. DYNAMICS AND STABILITY ANALYSIS 

In principle, the best-fit solutions found by the GAMP 
search should be dynamically stable. The chaos in the GM2 
solution can be explained by a relatively short integration 
time used in the algorithm (about 1000 Pc at 5 AU) which 
does not allow to eliminate mildly chaotic configurations. 
We assume here that the chaos manifestating through a lin- 
ear growth of MEGNO after a relatively long time (much 
longer than adopted in the GAMP) is regarded as a week one. 
This can be considered as a desired property of the method 
rather than its drawback. In fact, it can hardly be expected 
that multi-planet extrasolar systems found so far should be 
strictly quasi-periodic and regular, in terms of the Lyapunov 
exponent. Nevertheless, due to the fast convergence of the 
MEGNO and its great sensitivity to the chaotic evolution, 
the short integration times make it possible to eliminate the 
strongly chaotic configurations which very likely would lead 
to significant changes of the orbital elements. To verify the 
long-term stability of the best fit solutions found in this way, 
we should always analyze the structure of the phase space in 
their vicinity. Such analysis (using numerical or analytical 
tools) helps us to explore the dynamical environment of the 
nominal configuration, estimate whether the fits are robust to 
the errors and the found configurations do not change qualita- 
tively under small adjustments of the initial conditions. 

First, we have examinated the GMl fit by computing the 



stability map in the (fl(;,ec)-plane using the SM code. The rel- 
evant orbital parameters were varied and other elements were 
fixed at the values quoted in Table 1 . The results of this test 
are illustrated in Fig. |9] The integration time for this map is 
~ 10^ yr (~ 10"* orbital periods at 5 AU). The best fit solu- 
tion is marked by a big circle. The planetary collision line 
for fixed initial flb,eb is depicted as a smooth curve. The dy- 
namics over this line or in its proximity is very chaotic and 
the planets' eccentricities grow rapidly to 1, indicating col- 
lisions or ejections from the system. As we could expect, 
the stability map reveals an amazingly complex and beautiful 
structure of MMRs. These MMRs can be easily identified. It 
appears that the GMl fit is close to the border of the 1 1:2 and 
17:2 MMRs. In the previous paper on the HD 160691 system, 
we already suggested that this system may be very similar do 
the HD 12661 planetary system. It appears that, according 
to the current best-fit solutions, the dynamical environment of 
HD 160691 may be a close analogue of this system. Let us re- 
call that in Gozdziewski & Macieiewski (2003) we found that 
the HD 12661 may be located close to the 6: 1 MMR while the 
orbital fits of Butler et al. (2003) are related to the proximity 
of 1 1 :2 MMR ( Gozdziewski 2003 ). Due to the mentioned ten- 
dency that the estimated orbital period of the outermost body 
increases with a longer observational time span, the proxim- 
ity of the HD 160691 system to 6:1 MMR may also be quite 
possible. 

A similar stability test has been performed for the GM2 so- 
lution (Fig.llOk The upper panel is for the SM stability map 
in the (a ^,6^) -plane. The middle panel is for the maximum 
eccentricity of the planet c attained during the integration pe- 
riod. In both plots the best fit is marked with a big circle. For 
comparison, in the stability map (Fig. upper panel), the 
solution corresponding to NLl fit is marked with a smaller 
circle. Note that due to very similar orbital elements in the 
NLl and GM2 fits, depicting both of them at the same stabil- 
ity map is reasonable. Now it is clear why the NLl fit leads 
to a chaotic configuration. According to the identification of 
MMRs in this map, the GM2 fit is inside the libration island 
of the 5:1 MMR, while the NLl fit lies in a close vicinity or 
within the separatrix of this resonance. 

In the maxgc-map (Fig. [TO] the middle panel), we also 
marked (small circles) all fits gathered in the GAMP search 
within la confidence interval of the best fit solutions GMl, 
GM2. Such illustration of the best fit solutions is justified 
because, as we already mentioned, (flc,ec)-plane is represen- 
tative for the dynamical state of the system. The distribution 
of the best fits in this plane makes it possible to derive some 
interesting conclusions. In spite of the fact that the points rep- 
resenting the fits are spread over Oc £ [4,6) AU, the tendency 
of a decreasing with larger values of is clear Extrapo- 
lating this trend, the outer planet should have a small initial 
eccentricity for large Ac — 6-7 AU. The acceptable fit values 
within the la error of can be varied over 0-0.4. Similarly 
to the previously analyzed LNl and LN2 fits, the distribution 
of the osculating elements in the (ac,ec)-plane gives us an in- 
sight into realistic error estimates of these orbital elements. 

Obviously, all the relevant fits "avoid" the neighborhood of 
the collision zone. Close to and over the collision line, maxe^ 
grows to 1 which means that the system should be disrupted. 
The comparison of the stability maps for the GMl and GM2 
shows that the width of the MMRs significantly depends on 
the particular orbital parameters. Nevertheless, their positions 
do not change much. Let us recall, that the stability maps in 
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Fi g . l9l and 1 1 Olcorre spond to Uc altered by about 0.5 AU and 
by 2±0.1. 

The last panel shown in Fig. 1 101 refers to the maximum of 
the argument 9 = Ciic — C3b attained during the integration time. 
Librations of 9 correspond to the so called secular apsidal res- 
onance (SAR) when the planetary apsides remain aligned or 
anti-aligned. The behavior of 9 is important for the discussion 
of the dynamics of 2-planet configurations in the absence of 
strong MMRs. In this instance, the motion of the planetary 
system may be averaged over the mean longitudes. The secu- 
lar behavior of the averaged system can be described in terms 
of the s ecular octupole th eory of hierarchical planetary sys- 
tems bv lLee & Pealel ( l200 3b) and the recent more general sec- 
ular th eory of 2-planet systems by Michtchenko & Malhotra 
(|2004"). According to the results of these authors, the generic 
2-planet coplanar system which is far from low-order MMRs 
and collision zones exhibits a long-term stable behavior. This 
behavior can be characterized by the time evolution of 9. This 
angle can librate about 0° or 180°, circulate, or librate about 
0° in the large-eccentricity regime corresponding to the true 
secular resonance (for details see Michtchenko & Malhotra| 
1^04). When the eccentricities are moderate, the two first 
possibilities (modes) are alternative for the secular evolution 
of the p lanetary system. Under some conditions ( Lee & Peale 
|2003b), the SAR regime appears with a probability close to 1 . 

The octupole secular theory by Lee & Peale ( 2003b) is well 
suited for the HD 160691 system in the range of Oc > 4 AU 
and under the assumption that a given initial condition (IC) 
is distant from low-order MMRs. In this case, the best-fit 
solutions are characterized by a relatively small ratio a = 
db/flc — 0.3, as for hierarchical configurations. The bottom 
panel in Fig.fTOjillustrates the extent of the SAR in the range 
of flc e [4,6] AU. Besides the centers of the MMRs, the SAR 
with the anti-aligned apsides appears over an extended range 
of flc and for relatively small e^. Let us recall that this map has 
been derived numerically, simultaneously with calculating the 
SM stability map. It is very well reproduced by means of the 
octupole theory, see the map in the left panel in Fig.f^ which 
is obtained by th e numerical integration of the averaged equa- 
tions of motion (iLee & Pealell2003 b*a'). Another SAR map 
of the secular system, computed in the eccentricities plane, is 
shown in the right panel of Fig. In both of these maps, 
the best fit GMl is found in the extended zone of circulations 
of 9. Figure [T2l illustrates the evolution of 9(eb) and 9(ec) 
for the IC defined by GMl. These plots have been derived by 
the numerical integration of the full A^-body equations of mo- 
tion. The phase plots are computed for the constant values of 
the semi-major axes and the total angular momentum, J. Let 
us recall that the averaged system is governed by 1 degree of 
freedom Hamiltonian in which 9 and the eccentricity of one 
planet are the relevant canonical variables. To construct the 
phase plots, we fixed the initial 9 either at 0° or 180° and e\, 
was varied appropriately while the other initial eccentricity, 
gc, was derived from the J integral. The fact that the phase 
curves are not ideally smooth can be explained by the prox- 
imity to the two MMRs, 1 1:2 and 17:3. This confirms that the 
GMl fit corresponds to the configuration with 9 confined to 
the extended zone of circulation. The large amplitudes of ec- 
centricities shown in Fig.|8]can be also explained by the phase 
plots. 

The secular features of the GMl configuration remind us of 
another hierarchical system about HD 169830. The dynami- 
cal environments of their secular configurations are very sim- 
ilar. For instance, compare the discussed figures with the cor- 



responding ones from iGozdziewski & Konackil j2004h (their 
Fig. 14 and Fig. 15). The appearance of the SAR in extrasolar 
systems with well separated giant companions is not an ex- 
traordinary feature as sometimes claimed. In the light of the 
mentioned secular theories, in the regime of moderate eccen- 
tricities, the secular systems should be found in one of two 
distinct regimes of stable motion characterized by the circula- 
tion of librations of the critical argument 9. Nevertheless, one 
should be aware that the secular theories give us only an ap- 
proximate description of the long-term behavior. In the prox- 
imity to the strong MMRs, the stability should be investigated 
by an application of theories suitable for every individual res- 
onance or long-term direct integrations. 

5. CONCLUSIONS 

The new RV data of the HD 160691 s ystem published 
by the Anglo- Australian Planet Search team McCarth v et alJ 
i2004l) enable us to refine the orbital b est-fit solu tions. The 
preliminary observations published in iJones et a l. (2002^ 
suggested the presence of the second companion inducing a 
linear trend in the RV signal. In the light of the currently 
available data, it is very likely that the semi-major axis of the 
outer planet cannot be smaller than about 4 AU. This makes 
the system relatively separated. Curiously, the kinematic 2- 
Keplerian model of the RV signal constantly produces arti- 
facts e.g. the apparent proximity to the 2:1 MMR or large ec- 
centricities of the outer companion. These configurations lead 
to collisions between the planets in a very short time. In this 
work, we applied the "dynamical" fitting by forcing a require- 
ment that a real system should not be disrupted over a short 
time-span (thousands of years). Our analysis with the pro- 
posed GAMP algorithm reveals that the current data set still 
does not well constrain the semi-major axis of the outer com- 
panion. But the requirement of dynamical stability puts strong 
limits on the initial eccentricity of the outermost planet. The 
current set of observations is equally well modeled by a con- 
tinuum of statistically equivalent best-fit solution which are 
distributed over £ [4,6) AU and about [0 — 0.4], along a 
flat valley of (Xv)'''^ in the (ac,ec)-plane of osculating orbital 
parameters. Surprisingly, the initial longitudes of the coplanar 
configurations seem to be already well fixed. 

The dynamical analysis of the best-fit solutions reveals that 
the currently determined parameters of the outer planet are 
confined to the zone of low order MMRs with the inner planet, 
like 5:1, 11:2, 6:1, 13:2, 7:1. The distribution of the best so- 
lutions, within the formal la confidence interval of the best 
fit, reveals a clear correlation between initial values of and 
gc — the larger semi-major axis of the outer planet, the smaller 
its initial eccentricity. According to the empirically observed 
law, the increase of the observational window results in the 
increase of the determined orbital periods. This would mean 
that in the current best fit ~ 4.79 AU is still smaller than the 
real value. It is difficult to predict the real Oc but very likely it 
is in the range [5,6] AU. 

During the final preparation of this paper, we have 
learned about a discovery of a Neptunian-mass body in the 
HD 160691 system, close to the parent star dSantos et alJ 
^004). The semi-major axis of this small planet is about 
0.09 AU. Obviously, it is interesting to verify whether the 
presence of this body can affect the best-fit solutions and 
the dynamical behavior of the whole planetary system. The 
RV measurements published by the Anglo-Australian Planet 
Search team have the mean errors comparable to the RV semi- 
amplitude variations imposed by the new planet (4 m/s). For 
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the moment, the newest observations, as reported in the dis- 
covery paper, span a very short period of time (at present, 
about 80 d). Unfortunately, these measurements are not avail- 
able. They would be very helpful in improving our fits be- 
cause they cover the end-part of the longer period and should 
constrain the parameters of the outer companion. We tried to 
refine our fits using the data from a scanned figure published 
by the Swiss team on their WWW page^. The full data set, 
produced by a combination of the Anglo- Australian measure- 
ments and these synthetic RV points, has been firstly modeled 
using the A^-body signal. We used velocity offsets different for 
the two observatories. Still, we encountered similar problems 
as in the case of the Anglo-Australian data itself. The best 
A^-body fits for the 3 -planet system converge to Oc ~ 3.8 AU 
and gc > 0.6. According to the stability analysis, such fits 
locate the Jovian planets in a strongly chaotic zone. For the 3- 
planet best-fit solutions, we have (Xv)'''^ — 1-56 and the rms 
of ~ 3.6 m/s. Assuming that the innermost planet does not 
disturb the motion of the Jovian companions too much and ap- 
plying the GAMP for these two planets only (otherwise, the 
CPU time requirement is very significant), we derived solu- 
tions qualitatively similar to the GMl and GM2 fits (although 
the rms of these solutions exceeds 4.7 m/s). They tend to have 
a large ~ 0.5. Simultaneously, through the requirement of 
stabiHty, the GAMP places the best fits about Ac ~ 5.5 AU. 
Thus, in spite of a better time coverage, the curious inconsis- 
tency between the pure A^-body and the GAMP-derived solu- 
tions is still evident. These results can only be preliminary. 
A reliable modeling of the RV data is a difficult task without 
having access to the actual measurements. 

Although the elements of the innermost planet seem to be 
determined, due to the barely determined parameters of the 
outer planet, we think that the dynamical studies concerning 
the whole system do not make at present a great sense. As 
we have shown, the two Jovian planets can be coupled by a 
strong, low-order MMRs or can be well separated and rel- 
atively far from MMRs. These situations create completely 
different dynamical environments for the whole planetary sys- 
tem. 

In terms of the overall dynamical behavior, the /j Ara system 
appears to be similar to the HD 12661 and to the HD 169830. 



The Jovian planets in these systems are well separated. The 
variation of eccentricities spans a similar, relatively extended 
range from to 0.5-0.6. Their semi-major axes ratio corre- 
sponds to a proximity of low-order MMRs, like 6:1 or 11:2 
(for the HD 12661) or 7:1, 8:1 (similar to the HD 169830). 
Whether a real link exists between these systems is not clear 
at the moment. We still have to wait for new observations 
which eventually make it possible to derive significantly ac- 
curate determination of their parameters. 

Using the GAMP algorithm we performed a very prelimi- 
nary analysis of the RV data of the HD 154857 system. Sim- 
ilarity to the HD 160691 , these data reveal the presence of 
one Jovian planet and a linear trend, indicating the presence 
of the second, more distant body. Our A^-body fits make it 
possible to approximate the elements of this planet: its mass 
is about 0.5 Mj, semi-major axis about 4 AU and large ec- 
centricity ~ 0.78. Of course, these elements are barely con- 
strained. Nevertheless, the parameters of the inner planet are 
quite well fixed and we can akeady compute the stability map 
of this system. It is shown in Fig. ^] The relevant zone is 
filled up by strong MMRs, similarly to the HD 160691 sys- 
tem. By the same reasoning, as applied to the HD 160691 
system, the eccentricity of the outer planet should be strongly 
limited by the requirement of dynamical stability. 

Our work advocates a dynamical analysis of the RV data. 
Even if the observational window covers only partially the 
longest orbital period, we are able to derive a lot of charac- 
teristics of the discussed system, avoiding falsifying the or- 
bital elements by the purely kinematic approach. In fact, the 
analysis presented in this paper shows that the application of 
the GAMP-like algorithm is essential for finding reliable so- 
lutions which well reproduce the RV data and simultaneously 
lead to orbitally stable configurations of the planetary com- 
panions. 

6. ACKNOWLEDGMENTS 

This work is supported by the Polish Committee for Scien- 
tific Research, Grant No. 2P03D 001 22. K.G. is also sup- 
ported by the UMK grant 428-A. M.K. is also supported by 
NASA through grant NNG04GM62G. 



^ |http;//o~ 



qe.ch/'udry/planet/hdl60691.html| 



www.uni 



REFERENCES 



Adams, F. C. & Laughlin, G. 2003, Icarus, 163, 290 

Butler, R. R, Marcy, G. W.. Vogt, S. S., Fischer, D. A., Henry, G. W., 

Laughlin, G., & Wright, J. T. 2003, ApJ, 582, 455 
Butler, R. R, Tinney, C. G., Marcy, G. W., Jones, H. R. A., Penny, A. J., & 

Apps, K. 2001, ApJ, 555,410 
Charbonneau, R 1995, ApJS, 101, 309 

Cincotta, R M., Giordano, C. M., & Simo, C. 2003, Physica D, 182, 151 
Cincotta, R M. & Simo, C. 2000, A&AS, 147, 205 

Froeschle, C. & Lega, E. 1999, in NATO ASIC Proc. 522: The Dynamics of 
Small Bodies in the Solar System, A Major Key to Solar System Studies, 
463-+ 

Giordano, C. M. & Cincotta, R M. 2004, A&A, 423, 745 
Gozdziewski, K. & Konacki, M. 2004, ApJ, 610, 1093 
Gozdziewski, K., Konacki, M., & Maciejewski, A. J. 2003, ApJ, 594 
Gozdziewski, K. & Maciejewski, A. J. 2003, ApJ, 586, L153 
Gozdziewski, K. 2003, A&A, 398, 1 151 

Jones, H., Butler, R, Marcy, G., Tinney, C, Penny, A., C, M., & B., C. 2002a, 
MNRAS, astro-ph/0206216 



Jones, H. R. A., Raul Bufler, R., Marcy, G. W., Tinney, C. G.. Penny, A. J., 
McCarthy, C, & Carter, B. D. 2002b, MNRAS, 337, 1170 

Laskar, J. & Robutel, P. 1995, Celestial Mechanics and Dynamical 
Astronomy, 62, 193 

Laughlin, G. & Chambers, J. E. 2001, ApJ, 551, LI 09 

Lee, M. H. & Peale, S. J. 2003a, ApJ, 597, 644 

— . 2003b, ApJ, 592, 1201 

McCarthy, C. et al. 2004, ApJ 

Michtchenko, T. & Ferraz-Mello, S. 2001, ApJ, 122, 474 

Michtchenko, T. A. & Malhotra, R. 2004, Icarus, 168, 237 

Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. R 1992, 

Numerical Recipes in C. The Art of Scientific Computing (Cambridge 

Univ. Press) 
Santos, N. C. et al. 2004, A&A, 426, L19 

Stepinski, T. F, Malhotra, R., & Black, D. C. 2000, ApJ, 545, 1044 
Thommes, E. W. & Lissauer, J. J. 2003, ApJ, 597, 566 



7 



TABLE 1 

Initial conditions of the HD 160691 system. The first column is for 

THE BEST 2-KEPLER FIT EXPRESSED IN JACOBI ELEMENTS, ACCORDING TO THE 

model in Gozdziewski et al. 1 2003). Fit marked byNLI is derived by 
the self-consistent w-body fitting and expressed by osculating 

HELIOCENTRIC KEPLERIAN ELEMENTS AT THE MOMENT OF FIRST OBSERVATION 
{JD=2,450,915.291 1). THE NL2 SOLUTION IS DERIVED BY THE SAME 
APPROACH, WHEN THE MEASUREMENT ERRORS ARE SCALED BY THE STELLAR 
UTTER (~ 3 M/S). FITS DENOTED BY GMl AND GM2 ARE FORMALLY THE BEST 
FITS OBTAINED IN THE GAMP SEARCH. THE OSCULATING HELIOCENTRIC 
KEPLERIAN ELEMENTS ARE GIVEN AT THE DATE OF FIRST OBSERVATION. THE 
MASS OF THE PARENT STAR IS EQUAL TO 1.08 M ; . 



2K NLl NL2 GMl GM2 

Parameter be b c bcbcbc 

m2sin([Mj] 1.675 7.675 1.674(0.068) 2.707 (0.898) 1.683 5.689 1.670 2.981 1.677 2.398 

a[AU] 1.499 4.551 1.497 (0.007) 4.325 (0.560) 1.499 4.683 1.497 4.788 1.498 4.364 

P [d] 644.7 3368.7 

e 0.202 0.719 0.201 (0.029) 0.470(0.215) 0.203 0.585 0.204 0.340 0.201 0.411 

Q)[deg] 115.14 348.53 114.83(6.37) 341.14(16.23) 115.15 347.03 116.89 343.15 116.42 342.70 

Tp [JD-2,450,000] 2196.8 3681.4 

M[deg] 4.68 (4.54) 49.79(26.13) 3.67 66.50 1.45 65.2 2.75 46.42 

(Xiy/- 1.64378 1.64464 1.12871 1.64701 1.64707 

RMS [m/s] 4.15 4.15 4.12 4.17 4.15 

Vq [m/s] -37.2 -22.08 (6.12) -33.37 -22.08 -16.49 




IGOO L^GO 2D00 1500 [JD-2,4-50.000) [d] 



Fig. 1 . — Synthetic radial velocity curves for the HD 160691 system. The thin hne is for the best 2-Kepleria n solution (2K, see Tab le 1). The thick hne is for 
a stable, W-body fit GMl (see the text for explanation). Open circles are for the RV measurements published in IMcCarthv et all2004l) . 
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Fig. 2. — The best fits obtained by the Newtonian model for the RV data pubhshed in McCarthv et al. (2004). Small filled circles are for solutions with 
(Xv)'^^ within the formal 3a confidence interval of the NLl fit (see Table 1, (Xv)'^' < 1-75). Bigger, open circles are for (Xv)'^" < 1-667 corresponding to la 
confidence interval of the best fit NLl. Largest filled circles are for the fits of (Xv)'^" < 1-646, marginally larger, by about 1%, from the (Xv)'''^ — 1-644 of the 
best fit solution NLl. 




Fig. 3. — Orbital behavior in the best self-consistent, Newtonian fit NLl. Panel a) is for the eccentricities. Panel b) is for the argument of the SAR. Panel c) 
is for the semi-major axis of the outer companion. Panel d) is for the MEGNO, (7). The thin hne illustrates the Hnear fit = (Xniax/2)? + d. The derived 

Lyapunov time, Tl = \ /Xmax, is about 233 yr. 
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Fig. 4. — The best fits obtained by the Newtonian model of the RV data and their errors scaled by the stellar jitter with the mean dispersion about 3 m/s 
(see the text for explanation). Small filled circles are for solutions with (Xv)'^~ < 1-277, within the formal 3c confidence interval of the NL2 fit (see Table 1). 
Bigger, open circles are for (Zy)'^^ < 1.16 corresponding to la. Largest filled circles are for the fits of (x?)'^^ < 1130, marginally larger, by about 1%, from 
the ixiy^^ ^ 1-129 ofNL2. 
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Fig. 5. — A comparison of the disfiibution of the best fits obtained by the A'-body fitting in the (ac,ec)-plane. The left panel is for the original data set, the 
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Fig. 6. — The best fits obtained by the GAMP. Open circles are for (Xv)''^^ < 1-668, i.e., within formal la confidence interval of the best fit solution having 
= 1-647. 




Fig. 7. — Orbital behavior in the GM2 fit. Panel a) is for the eccentricities. Panel b) is for the argument of the SAR. Panel c) is for the critical argument of the 
5:1 MMR. Panel d) is for the MEGNO, (Y). The thin line illustrates the linear fit (Y){t) = (Vax/2)f + d. The derived Lyapunov time, 7i = 1/Vax, is about 
35,000 yr. 
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Fig. 8. — Orbital behavior in the GMl fit. Panel a) is for the eccentricities. 
Y{t). Panel d) is for the MEGNO, (F). 
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Fig. 9. — The stability map in the (ac.fcj-plane, in terms of the spectral number, logSN, (see the text for an explanation) for the best fit solution GMl (see 
Table 1). Color code is for logS'A': black means quasi-periodic, regular configurations, and yellow strongly chaotic systems. Labels are for the identification of 
the low-order MMRs. A big circle marks the parameters of the GMl fit. 
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Fig. 10. — The upper plot is for the stability map in the (ac,ec)-plane, in term.s of the spectral number, logSN, (see the text for an explanation) for the best 
fit solution GM2 (see Table 1). Color code is for logSW: black means quasi-periodic, regular configurations and yellow strongly chaotic systems. A big circle 
marks the parameters of the GM2 fit. For comparison, the elements of NLl fit are depicted as a smaller circle. The middle plot is for maximum eccentricity of 
the outer planet attained during the integration period (about lO' yr). Circles in this plot mark parameters of the fits obtained in the GAMP search within lo 
confidence interval of the GMl fit. The lower plot is for the maximum of 0. The resolution of the maps is is 600 x 100 data points. 
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Fig. 1 1 . — The semi-amplitude of the critical ai'gument of SAR, 9, about the libration center 180°, derived by the application of the secular oc tupo le theoiy by 
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Fig. 12. — Phase diagram of the GMl fit obtained by numerical integration of the full, W-body model for constant values of the integrals of the total angular 
momentum and total energy. The nominal configuration is marked by a filled circle in both panels. See the text for explanation. 
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Fig. 13. — The stability map in the (acj^c) -plane, in terms of the spectral number, SN, for the HD 154857 planetary system. Color code is for SN: black means 
quasi-periodic, regular configurations and yellow (light gray) strongly chaotic systems. Labels are for the identification of the low-order MMRs. Resolution of 
the plot is 640 x 100 data points. 



